Texture Analysis

GLCM : Gray-level Co-occurrence Matrices

The Gray-Level Co-occurrence Matrix (GLCM) counts the distribution of gray values ​​between pixels in a grayscale image to distinguish different textures.

Gray Level Co-occurrence matrix (GLCM) uses adjacency concept in images. The basic idea is that it looks for pairs of adjacent pixel values that occur in an image and keeps recording it over the entire image. Below figure explains how a GLCM is constructed.

s you can see from the above image, gray-level pixel value 1 and 2 occurs twice in the image and hence GLCM records it as two. But pixel value 1 and 3 occurs only once in the image and thus GLCM records it as one. Of course, I have assumed the adjacency calculation only from left-to-right. Actually, there are four types of adjacency and hence four GLCM matrices are constructed for a single image. Four types of adjacency are as follows.

• Left-to-Right
• Top-to-Bottom
• Top Left-to-Bottom Right
• Top Right-to-Bottom Left

The degree co-occurrence matrix can be defined as a gray scale as i i The gray value of the pixel point and another pixel point corresponding to it is j j The probability. Then all estimated values ​​can be expressed in the form of a matrix, which is called a gray-level co-occurrence matrix. Such as: according to any point in the image (x,y) ( x , y ) The gray value of and its corresponding point (x+dx,y+dy) ( x + d x , y + d y ) The gray value of can get a gray value combination (g1,g2) ( g 1 , g 2 ) . Count the probability matrix of each combination of gray values ​​in the whole image P P It is the gray level co-occurrence matrix.

Due to the large dimensionality of the gray level co-occurrence matrix, it is generally not directly used as a feature to distinguish textures, but based on some statistics constructed by it as texture classification features. 14 kinds of statistics calculated based on the gray level co-occurrence matrix have been proposed: energy, entropy, contrast, uniformity, correlation, variance, sum average, sum variance, sum entropy, difference variance, difference average, difference entropy, related information Measure and maximum correlation coefficient.

GLCM Features

• Contrast (CON): The contrast measures the spatial frequency of an image and is a different moment of GLCM. It is the difference between the highest and the lowest values of the adjacent set of pixels. The contrast texture measures the local variations present in the image. An image with the low contrast presents GLCM concentration term around the principal diagonal and features low spatial frequencies.

• Homogeneity (HOM): This statistical measure is also called inverse difference moment. It measures the homogeneity in the image where it assumes larger values for smaller differences in grey tone within-pair elements. Homogeneity is more sensitive to the presence of near diagonal elements in the GLCM. The value of homogeneity is maximum when elements in the image are the same. GLCM contrast and homogeneity are strongly but inversely correlated, which means homogeneity decreases when contrast increases while energy is kept constant.

• Dissimilarity (DIS): Dissimilarity is a linear measure of local variations in an image.

• Angular second Moment (ASM): It measures textural uniformity i.e. repetitions in pixel pair. It detects the disorders in textures of the images. The maximum value achieved by the angular second moment is one. Higher values occur when the gray level distribution has a constant periodic form.

• Energy (EG): Energy is computed as the square root of an angular second moment. When the window is orderly, energy has higher values.

• Correlation (CORR): It is a measure of linear dependencies between the gray tone of the image.

Let's download a Sentinel 2 image using the Google Earth Engine for this example:

!pip install rasterio

Collecting rasterio
  Downloading rasterio-1.3.11-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (14 kB)
Collecting affine (from rasterio)
  Downloading affine-2.4.0-py3-none-any.whl.metadata (4.0 kB)
Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio) (24.2.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio) (2024.8.30)
Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from rasterio) (8.1.7)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from rasterio) (0.7.2)
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.26.4)
Collecting snuggs>=1.4.1 (from rasterio)
  Downloading snuggs-1.4.7-py3-none-any.whl.metadata (3.4 kB)
Requirement already satisfied: click-plugins in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.1.1)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio) (71.0.4)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio) (3.1.4)
Downloading rasterio-1.3.11-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (21.7 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 21.7/21.7 MB 25.7 MB/s eta 0:00:00
Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)
Downloading affine-2.4.0-py3-none-any.whl (15 kB)
Installing collected packages: snuggs, affine, rasterio
Successfully installed affine-2.4.0 rasterio-1.3.11 snuggs-1.4.7

import ee
ee.Authenticate()
ee.Initialize(project='my-project-1527255156007')

from google.colab import drive
drive.mount('/content/drive')

Mounted at /content/drive

import folium
from folium import plugins
from IPython.display import Image
import rasterio
import numpy as np
from matplotlib import pyplot as plt
import cv2

AOI = ee.Geometry.Polygon(
        [[[2.104356850541711, 41.44948070383667],
          [2.104356850541711, 41.3574829950726],
          [2.250869835771203, 41.3574829950726],
          [2.250869835771203, 41.44948070383667]]])

startDateviz = ee.Date.fromYMD(2021,1,1);
endDateviz = ee.Date.fromYMD(2021,1,31);
collectionviz = ee.ImageCollection("COPERNICUS/S2_HARMONIZED").filterDate(startDateviz,endDateviz).filterBounds(AOI).filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 10);

S2 = collectionviz.median().clip(AOI).divide(10000)

image = S2.select(['B2','B3','B4','B5','B8']).reproject('EPSG:4326', None, 20)

task = ee.batch.Export.image.toDrive(image=image,
                                         crs='EPSG:4326',
                                         scale=20,
                                         fileFormat='GeoTIFF',
                                         description='Barcelona' ,
                                         maxPixels=1e13,
                                         region= AOI)

task.start()

Now let's open the image stored in Google Drive:

path_barcelona = '/content/drive/MyDrive/Datasets_CV/Barcelona.tif'

with rasterio.open(path_barcelona) as src:
    im = src.read()

im = im.transpose([1,2,0])

img = im*2*255

img = img.astype('uint8')

red = img[:,:,2]
green = img[:,:,1]
blue = img[:,:,0]

And plot the RGB image:

rgb = np.dstack((red, green, blue))
plt.figure(figsize=[12,12])
plt.imshow(rgb)

<matplotlib.image.AxesImage at 0x7da66e0c7430>

Let's convert the image to HSV and use the intensity represented by the Value band:

bgr = np.dstack((blue, green, red))

img_hsv = cv2.cvtColor(bgr, cv2.COLOR_BGR2HSV)

First plot the converted image to HSV:

plt.figure(figsize=[12,12])
plt.imshow(img_hsv)

<matplotlib.image.AxesImage at 0x7da66d002440>

And then only Value band on 'inferno' colormap:

img_value = img_hsv[:,:,2]

plt.figure(figsize=[12,12])
plt.imshow(img_value, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da6af091480>

GLCM from skimage

In this first part we will select regions in the image that represent two different types of LULC. We will select 3 regions of 16 by 16 pixels for the Urban class and 3 for the Water class.

from skimage.feature import graycomatrix, graycoprops

PATCH_SIZE = 16
water_locations = [(454, 600), (350, 700), (290, 750)]
water_patches = []
for loc in water_locations:
    water_patches.append(img_value[loc[0]:loc[0] + PATCH_SIZE,
                               loc[1]:loc[1] + PATCH_SIZE])

urban_locations = [(400, 200), (300, 300),(120, 510)]
urban_patches = []
for loc in urban_locations:
    urban_patches.append(img_value[loc[0]:loc[0] + PATCH_SIZE,
                             loc[1]:loc[1] + PATCH_SIZE])

Let's plot the location of these regions on our image:

fig = plt.figure(figsize=(10, 10))

# display original image with locations of patches
ax = fig.add_subplot(1, 1, 1)
ax.imshow(img_value, cmap='inferno',
          vmin=0, vmax=255)
for (y, x) in water_locations:
    ax.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'gs')
for (y, x) in urban_locations:
    ax.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'bs')

Now let's take a closer look at how these images are represented in the Intensity values:

fig = plt.figure(figsize=(10, 10))
for i, patch in enumerate(water_patches):
    ax = fig.add_subplot(3, len(water_patches), len(water_patches)*1 + i + 1)
    ax.imshow(patch, cmap=plt.cm.gray,
              vmin=0, vmax=255)
    ax.set_xlabel('Water %d' % (i + 1))

for i, patch in enumerate(urban_patches):
    ax = fig.add_subplot(3, len(urban_patches), len(urban_patches)*2 + i + 1)
    ax.imshow(patch, cmap=plt.cm.gray,
              vmin=0, vmax=255)
    ax.set_xlabel('Urban %d' % (i + 1))
plt.show()

We saw that for class Water the pixel intensity is homogeneous throughout the selected region, unlike the urban area, where different intensity values are found in the same region.

Knowing this difference in the texture of the two classes, let's apply the GLCM and get the values of only two characteristics: Dissimilarity and Correlation. In the skimage function we will just use a value on pixel distance and an angle.

xs = []
ys = []
for patch in (water_patches + urban_patches):
    glcm = graycomatrix(patch, distances=[5], angles=[0], levels=256,
                        symmetric=True, normed=True)
    xs.append(graycoprops(glcm, 'dissimilarity')[0, 0])
    ys.append(graycoprops(glcm, 'correlation')[0, 0])

So we can plot a graph relating the two properties of each region:

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.plot(xs[:len(water_patches)], ys[:len(water_patches)], 'go',
        label='Water')
ax.plot(xs[len(urban_patches):], ys[len(urban_patches):], 'bo',
        label='Urban')
ax.set_xlabel('GLCM Dissimilarity')
ax.set_ylabel('GLCM Correlation')
ax.legend()
plt.show()

Let's now increase the amount of angles and displacement distances, thus obtaining a matrix with the values referring to the amount of angles and distances added:

glcm_urban = graycomatrix(urban_patches[0], distances=[1,2], angles= [0, np.pi/4, np.pi/2, 3*np.pi/4], levels=256,symmetric=True, normed=True)

So we can get all the GLCM characteristics of the selected region:

contrast = graycoprops(glcm_urban, 'contrast')
dissimilarity = graycoprops(glcm_urban, 'dissimilarity')
homogeneity = graycoprops(glcm_urban, 'homogeneity')
energy = graycoprops(glcm_urban, 'energy')
correlation = graycoprops(glcm_urban, 'correlation')
asm = graycoprops(glcm_urban, 'ASM')

print('contrast:', contrast)

contrast: [[219.075      719.72       601.98333333 649.05777778]
 [428.58035714 719.72       816.16964286 649.05777778]]

print('dissimilarity:', dissimilarity)

dissimilarity: [[ 8.98333333 17.37333333 15.725      16.89777778]
 [12.30357143 17.37333333 17.19642857 16.89777778]]

print('homogeneity:', homogeneity)

homogeneity: [[0.14944188 0.08130337 0.10498005 0.09038129]
 [0.14221037 0.08130337 0.10876507 0.09038129]]

print('energy:', energy)

energy: [[0.05818433 0.05803362 0.05667279 0.05708992]
 [0.05734479 0.05803362 0.05803571 0.05708992]]

print('correlation:', correlation)

correlation: [[0.79179501 0.34858486 0.42746833 0.4128382 ]
 [0.59604784 0.34858486 0.13545407 0.4128382 ]]

print('asm:', asm)

asm: [[0.00338542 0.0033679  0.00321181 0.00325926]
 [0.00328842 0.0033679  0.00336814 0.00325926]]

Using fast_glcm:

The implementation of skimage results in a single value for each selected region. To get an image representing each GLCM property, we're going to use an implementation that runs through the entire original image extracting texture information locally. This way we will use the fast_glcm (link on reference) implementation to find the GLCM properties of the complete image:

mi, ma = 0, 255
nbit=8
ks = 5
h,w = img_value.shape

    # digitize
bins = np.linspace(mi, ma+1, nbit+1)
gl1 = np.digitize(img_value, bins) - 1
gl2 = np.append(gl1[:,1:], gl1[:,-1:], axis=1)

    # make glcm
glcm = np.zeros((nbit, nbit, h, w), dtype=np.uint8)
for i in range(nbit):
    for j in range(nbit):
        mask = ((gl1==i) & (gl2==j))
        glcm[i,j, mask] = 1

kernel = np.ones((ks, ks), dtype=np.uint8)
for i in range(nbit):
    for j in range(nbit):
        glcm[i,j] = cv2.filter2D(glcm[i,j], -1, kernel)

glcm = glcm.astype(np.float32)

Now we will plot the resulting image for each property:

Contrast

cont = np.zeros((h,w), dtype=np.float32)
for i in range(nbit):
    for j in range(nbit):
        cont += glcm[i,j] * (i-j)**2

plt.figure(figsize=[12,12])
plt.imshow(cont, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da66147abc0>

Dissimilarity

diss = np.zeros((h,w), dtype=np.float32)
for i in range(nbit):
    for j in range(nbit):
        diss += glcm[i,j] * np.abs(i-j)

plt.figure(figsize=[12,12])
plt.imshow(diss, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da66122b160>

Homogeneity

homo = np.zeros((h,w), dtype=np.float32)
for i in range(nbit):
    for j in range(nbit):
        homo += glcm[i,j] / (1.+(i-j)**2)

plt.figure(figsize=[12,12])
plt.imshow(homo, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da6610ba2c0>

ASM

asm = np.zeros((h,w), dtype=np.float32)
for i in range(nbit):
    for j in range(nbit):
        asm  += glcm[i,j]**2

plt.figure(figsize=[12,12])
plt.imshow(asm, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da660f40be0>

Energy

ene = np.sqrt(asm)

plt.figure(figsize=[12,12])
plt.imshow(ene, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da660faf640>

Correlation

pnorm = glcm / np.sum(glcm, axis=(0,1)) + 1./ks**2
ent  = np.sum(-pnorm * np.log(pnorm), axis=(0,1))

plt.figure(figsize=[12,12])
plt.imshow(ent, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da6610277c0>

Local Binary Patterns

There are lots of different types of texture descriptors are used to extract features of an image. Local Binary Pattern, also known as LBP, is a simple and grey-scale invariant texture descriptor measure for classification. In LBP, a binary code is generated at each pixel by thresholding it’s neighbourhood pixels to either 0 or 1 based on the value of the centre pixel.

The rule for finding LBP of an image is as follows:

Set a pixel value as center pixel. Collect its neighbourhood pixels (Here I am taking a 3 x 3 matrix so; total number of neighbourhood pixel is 8) Threshold it’s neighbourhood pixel value to 1 if its value is greater than or equal to centre pixel value otherwise threshold it to 0. After thresholding, collect all threshold values from neighbourhood either clockwise or anti-clockwise. The collection will give you an 8-digit binary code. Convert the binary code into decimal. Replace the center pixel value with resulted decimal and do the same process for all pixel values present in image. Let’s take an example to understand it properly.

Let’s take a pixel value from the above output to find its binary pattern from its local neighbourhood. So, I am taking a value ‘149’ (present at 15th row and 19nd column) and its 8 neighbourhood pixels to form a 3 x 3 matrix.

Collect the thresholding values either clockwise or anti-clockwise. Here, I am collecting them clockwise from top-left. So, after collecting, the binary value will be as follows:

Then, convert the binary code into decimal and place it at center of matrix.

1 x 27 + 1 x 26 + 1 x 25 + 0 x 24 + 0 x 23 + 0 x 22 + 0 x 21 +1 x 20 = 128 + 64 + 32 + 0 + 0 + 0 + 0 + 1 = 225

Now, the resulted matrix will look like,

Let's use the skiimage implementation and apply it to our test image:

from skimage.feature import local_binary_pattern

Let's define the two parameters: the radius and the number of neighbors around the pixel:

radius = 1
n_points = 8 * radius

patterns = local_binary_pattern(img_value, n_points, radius, 'uniform')

We can then plot the result of applying the function to the image:

plt.figure(figsize=[12,12])
plt.imshow(patterns, cmap = 'inferno' )

<matplotlib.image.AxesImage at 0x7da660eaf4f0>

The image can be converted into a histogram, offering a better representation of its textures for classification using some Machine Learning algorithm:

ax = plt.hist(patterns.flatten(), bins = 256)
plt.show()

References:

https://www.programmersought.com/article/96625951757/

https://gogul.dev/software/texture-recognition

https://scikit-image.org/docs/dev/auto_examples/features_detection/plot_glcm.html

https://scikit-image.org/docs/0.7.0/api/skimage.feature.texture.html

https://peerj.com/articles/cs-536/

https://www.geeksforgeeks.org/create-local-binary-pattern-of-an-image-using-opencv-python/

https://github.com/tzm030329/GLCM

https://www.linkedin.com/pulse/image-classification-opencv-scikit-learn-raj-ramanujam/

https://www.geeksforgeeks.org/create-local-binary-pattern-of-an-image-using-opencv-python/

https://fairyonice.github.io/implement-lbp-from%20scratch.html